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ABSTRACT 0 6 / ' - - 3 

We present the properties of accretion disk corona (ADC) models in which the radiation field, the 
temperature, and the total opacity of the corona are determined self-consistently. We use a nonlinear 
Monte Carlo code to perform the calculations. As an example, we discuss models in which the corona is 
situated above and below a cold accretion disk with a plane-parallel (slab) geometry, similar to the 
model of Haardt & Maraschi. By Comptonizing the soft radiation emitted by the accretion disk, the 
corona is responsible for producing the high-energy component of the escaping radiation. Our models 
include the reprocessing of radiation in the accretion disk. Here the photons either are Compton- 
reflected or photoabsorbed, giving rise to fluorescent line emission and thermal emission. The self- 
consistent coronal temperature is determined by balancing heating (due to viscous energy dissipation) 
with Compton cooling, determined using the fully relativistic, angle-dependent cross sections. The total 
opacity is found by balancing pair productions with annihilations. We find that, for a disk temperature 
kT BB < 200 eV, these coronae are unable to have a self-consistent temperature higher than ~ 140 keV if 
the total optical depth is >0.2, regardless of the compactness parameter of the corona and the seed 
opacity. This limitation corresponds to the angle-averaged spectrum of escaping radiation having a 
photon index >1.8 within the 5-30 keV band. Finally, all models that have reprocessing features also 
predict a large thermal excess at lower energies. These constraints make explaining the X-ray spectra of 
persistent black hole candidates with ADC models very problematic. 

Subject headings: accretion, accretion disks — radiation mechanisms: nonthermal — radiative transfer 
stars: coronae 
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1. INTRODUCTION 

The high-energy spectra from Galactic black hole candi- 
dates (BHCs) and radio-quiet active galactic nuclei (AGNs) 
are similar. Both can be described roughly by a power law, 
with a photon index of ~ 1.5 2.0 that extends up to — 100 
keV, modified by an exponential cutoff (Grebenev et al. 
1993; Gilfanov et al. 1993; Maisack et al. 1993; Wilms et al. 
1996, and references therein). The spectra of many compact 
objects also show features due to reprocessing of X-rays 
with moderately cold matter, including a fluorescence line 
due to iron and a Compton reflection hump (Guilbert & 
Rees 1988; Pounds et al. 1990; Nandra et al. 1991; Fabian 
1994, and references therein). 

Comptonization of radiation by a semirelativistic plasma, 
in the form of a corona associated with an accretion disk, 
may be the principal mechanism that produces the high- 
energy portion of these spectra. However, even within the 
context of these models, the geometry and physical proper- 
ties of the corona are still being debated. One popular 
model involves an optically thick, cold (kT < 1 keV) accre- 
tion disk and a hot accretion disk corona (ADC), often 
having a plane-parallel (slab) configuration (Galeev, 
Rosner, & Vaiana 1979; Sunyaev & Titarchuk 1980; Field 
& Rogers 1993; Haardt & Maraschi 1991, 1993, hereafter 
HM91, HM93, respectively; Haardt, Maraschi, & Ghisellini 
1994, 1997; Nakamura & Osaki 1993; Hua & Titarchuk 
1995; Stern et al. 1995a; Poutanen & Svensson 1996, and 
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references therein). One nice feature of the ADC models is 
that the accretion disk serves a double role in being the 
source of the seed photons that are Comptonized by the 
corona and being the medium responsible for the reprocess- 
ing and reflection features. It is commonly assumed that a 
substantial fraction of the gravitational accretion energy is 
dissipated directly into the corona, although the physical 
processes by which the corona is heated are still unknown. 

Recent theoretical work provides a heuristic framework 
for understanding why accretion disk coronae should form. 
Balbus & Hawley (1991; see also Hawley, Gammie, & 
Balbus 1994 and references therein) have found a powerful 
MHD instability that, in its nonlinear development, tends 
to amplify the turbulent magnetic pressure to a level on the 
order of the gas pressure. Buoyancy of the magnetic flux 
tubes will tend to make them rise out of the disk’s interior 
(Galeev et al. 1979; Tout & Pringle 1992), into a region 
where they dominate the ambient gas pressure. As is 
believed to occur in the solar corona, this type of magnetic 
field evolution gives rise to a field structure in which large- 
scale magnetic reconnection is inevitable, releasing energy 
that can heat the corona (Hawley et al. 1997; Balbus, 
Hawley, & Stone 1996). More recent work has shown that 
the magnetic field energy density is increased in low-density 
regions (i.e., the outer disk atmosphere), possibly allowing 
for direct coronal heating (Stone et al. 1996). 

Alternate accretion disk models also have been able to 
explain the high-energy spectra of X-ray binaries. One 
common alternative model involves an optically thin, hot 
accretion disk (Shapiro & Lightman 1976; Kusunose & 
Mineshige 1995; White & Lightman 1990; Luo & Liang 
1994; Narayan & Yi 1995 and references therein; Chen 
1995 and references therein). However, these models, unlike 
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the ADC models, must include an additional reprocessing 
medium (e.g., a stellar wind or an outer accretion disk) in 
order to explain the observed reflection features. 

Because of the coupling between the radiation field and 
the coronal plasma, solving the nonlinear radiative transfer 
problem for accretion disk coronae is very complicated. 
Through Comptonization, the properties of the radiation 
field depend on the corona's temperature and opacity. Since 
Compton cooling is the dominant cooling process, the tem- 
perature of the corona, in turn, depends on the radiation 
field. In addition, the corona and radiation field are coupled 
to each other through the processes of pair production and 
pair annihilation. Reprocessing of radiation within the acc- 
retion disk, where the radiation is either Compton-reflected 
back into the corona or photoabsorbed, also complicates 
the problem. 

For accretion disk corona models, the high-energy spec- 
trum of escaping radiation depends primarily on the 
geometry of the system, the Thomson optical depth of the 
corona, t t , and the temperature of the corona, T c . For a 
given geometry, the last two parameters typically are 
allowed to vary independently until a good fit to the data is 
achieved. However, an ADC system with a given T c and t x 
is not necessarily a self-consistent model. For many geome- 
tries, reprocessing of radiation from the corona within the 
cold accretion disk allows the accretion disk to have a large 
flux of thermal radiation even if most of the gravitational 
energy is dissipated within the corona (HM91). As we will 
show for the specific case of slab geometries, a high flux of 
soft photons results in high Compton cooling rates, and 
consequently, models with modest optical depths (t x > 0.2) 
can only have temperatures T c < 120 keV. Therefore, pre- 
vious ADC models that yield acceptable fits to the data of 
BHCs (Haardt et al. 1993; Titarchuk 1994) are not physi- 
cally realizable. For example, the spectrum of Cyg X-l has 
recently been described by an ADC model, with a slab 
geometry, having a temperature of T c ~ 1 50 keV and an 
optical depth of i T - 0.3 (Haardt et al. 1993). As we discuss 
below, this is not a self-consistent combination. 

Haardt & Maraschi (HM91, HM93; Haardt 1993) devel- 
oped the first thermal Comptonization model in which the 
coupling between the radiation field and the coronal plasma 
is taken into account. Reprocessing of radiation within the 
cold accretion disk also is included in their model, and a 
self-consistent temperature of the corona is determined by 
equating Compton cooling with viscous heating. As we will 
discuss in § 2, this model uses a hybrid of Monte Carlo (MC) 
and analytic techniques. However, several of the approx- 
imations used in calculating the radiation field limit the 
accuracy and validity of these calculations. 

In this paper, we present an ADC model that uses a 
nonlinear Monte Carlo (NLMC) code adapted from Stern’s 
"‘large particle” Monte Carlo technique (Stern 1985; Stern 
et al. 1995a, and references therein). Although we only con- 
sider the slab geometry in this paper, more complicated 
geometries can be modeled and will be the focus of future 
papers. Rather than propagating photons through a 
“background” medium one at a time, as is the case with 
linear Monte Carlo methods (Pozdnyakov, Sobol, & 
Sunyaev 1983), all types of particles are propagated in 
parallel (for this paper, the types of particles included are 
photons, electrons, and positrons). In addition, the ener- 
getics of the system, including the Compton cooling rate, 
the pair production rate, and the electron/positron annihi- 


lation rate, is calculated numerically. Thus the pair opacity, 
the temperature of the corona, and the radiation field can be 
computed self-consistently. These quantities are computed 
locally , giving rise to models with nonuniform temperature 
and/or density distributions (see discussion in § 5). Fully 
relativistic and angle-dependent cross sections are used, and 
there is no need to make any radiative transfer approx- 
imations. The pair production rates and Compton cooling 
rates are by computed using the self-consistent, angle- 
dependent radiation field rather than analytic approx- 
imations. 

The remainder of the paper is organized as follows : In § 2, 
we discuss the limitations of the computational methods 
used in previous ADC models. In § 3, we summarize the 
“large particle” Monte Carlo method and describe our 
modifications that make it more applicable to thermal acc- 
retion disk corona models. In § 4, we describe our ADC 
model and note the differences between it and previous 
models. In § 5, we present the results of our simulations, and 
in § 6 we summarize our results and give our plans for future 
work. 

2. PREVIOUS ADC MODELS AND COMPUTATIONAL 
TECHNIQUES 

To date, most thermal Comptonization models have used 
either the analytic methods pioneered by Lightman & 
Rybicki (1980) and Sunyaev & Titarchuk (1980) (and 
expanded upon by Sunyaev & Titarchuk 1985 and Hua & 
Titarchuk 1995), numerical methods of solving the kinetic 
equations (Zdziarski 1985; Svensson 1987; Lightman & 
Zdziarski 1987; Coppi 1992; Poutanen & Svensson 1996, 
and references therein), or “linear” Monte Carlo (MC) 
methods (Pozdnyakov et al. 1983; Gorecki & Wilczewski 
1984). Since the properties of the corona have to be 
assumed before the radiation field is computed with an 
analytic method, it is uncertain whether the properties of 
the corona are self-consistent with the radiation field. In 
addition, analytic solutions can be found for only a small 
number of simple coronal geometries (most importantly, 
spherical or plane-parallel geometry), and it does not 
appear possible to perform the computations using the fully 
relativistic, angle-dependent cross sections (Hua & Titar- 
chuk 1995). Also, no analytic method has been presented 
that models the reprocessing of radiation in the accretion 
disk while taking into account the anisotropy of the radi- 
ation field within the corona. Thus, if reprocessing is impor- 
tant in the source, “reflection features” have to be 
computed by using a different computational method and 
added to the escaping radiation field. With this approach, 
subsequent interactions between the reflected component of 
the radiation field and the ADC cannot properly be taken 
into account. 

Linear MC methods, in which one photon at a time is 
propagated through a fixed medium, have to specify the 
properties of the corona before simulating the radiation 
field ; as is the case with the analytic models, the specified 
properties of the corona often are not self-consistent. In 
addition, since linear MC methods follow one photon at a 
time, simulating photon-photon pair productions and anni- 
hilations must be done analytically, and pair opacity must 
be approximated by interfacing these semianalytic calcu- 
lations with the MC computations (HM93). 

As discussed in § 1, Haardt & Maraschi (HM91, HM93) 
developed an improved thermal Comptonization model for 
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plane-parallel accretion disk coronae by taking the coup- 
ling between the accretion disk and the corona into 
account. The temperature of the corona is calculated by 
balancing the total Compton cooling rate with the assumed 
heating rate (which is determined by the coronal com- 
pactness, a free parameter). All the gravitational accretion 
energy is dissipated uniformly into the corona. Therefore all 
the energy emitted by the cold accretion disk is supplied by 
irradiation from the corona. The total opacity is calculated 
by balancing photon-photon pair productions with annihi- 
lations. 

The models by HM91 and HM93 are the first models in 
which the coronal temperature and opacity are not inde- 
pendent free parameters. However, because of the analytic 
approximations, these models have several shortcomings: 

(1) Since HM93 use the Thomson cross section in the elec- 
tron’s rest frame rather than the correct Klein-Nishina cross 
section, a Wien cutoff, rather than a self-consistent spectral 
shape, is used to approximate the spectrum for energies 
e > kT c , where T c is the temperature of the corona. This 
approximation leads to an underestimation of the pair pro- 
duction rates since the high-energy tail of the self-consistent 
spectrum is harder than a Wien cutoff (Stern et al. 1995b). 

(2) Additional uncertainties arise because the Compton 
cooling rate is calculated by approximating the spectrum of 
radiation as a power law with a Wien cutoff. (3) HM93 
neglect multiple reflections, line features, and subsequent 
interactions with the corona by reflected photons. Very 
recently, Haardt et al. (1997) have modified the HM93 
model such that the high-energy rollover is calculated more 
accurately by using a fully relativistic kernel for isotropic 
unpolarized radiation (see Haardt et al. 1997 for a detailed 
discussion). 

Using an iterative scattering method, Poutanen & 
Svensson (1996) solved the radiative transfer problem for 
ADC models, in which Compton scattering, photon-photon 
pair production, pair annihilation, bremsstrahlung, and 
double Compton scattering are taken into account. In addi- 
tion, the relativistic cross sections are used, the spectrum of 
escaping radiation can be determined for any direction, and 
the reprocessing of coronal radiation in the accretion disk is 
taken into account by using a reflection matrix. The main 
limitation of this method is that it is one-dimensional, and 
complicated geometries therefore cannot be considered. 
Another limitation is that solutions can be found only for 
small optical depths (t t < 1, depending on the coronal tem- 
perature; see Poutanen & Svensson 1996 for more details). 

3. THE NONLINEAR MONTE CARLO CODE 
3.1. Overview 

To study the nonlinear problem of radiative transfer in 
an ADC, we have modified the "large particle” Monte 
Carlo (LPMC) method, developed by Stern (1985, 1988) 
and Stern et al. (1995a, 1995b). Unlike most Monte Carlo 
codes used for accretion disk coronae, in which photons are 
propagated through a "background” medium one at a 
time, the LPMC code simulates all particles (photons, elec- 
trons, positrons) in parallel (i.e., simultaneously). The 
advantage of propagating all particles in parallel is that 
there is no "background” medium; photons can interact 
with electrons, and like particles can interact with each 
other. Because the coupling between the radiation field and 
the plasma is taken into account by directly simulating all 


particles in parallel, the nonlinear radiation transfer 
problem is solved self-consistently. 

One challenge of using Monte Carlo numerical tech- 
niques for simulating X-ray spectra is to obtain good sta- 
tistics with models requiring only a modest amount of CPU 
time. Since we are simulating high-resolution spectra that 
span many decades in energy, it is not possible to simulate 
the particles by just interpreting any particle within the 
code as a particle in the real physical system. If this 
approach were used, a great deal of CPU time would be 
wasted by tracking the low-energy portion of the photon 
spectrum. To deal with this problem. Stern (1985) developed 
a numerical technique using “large particles” (LPs). Each 
LP represents a group of identical particles having the same 
particle type (photon, electron, or positron), energy, veloc- 
ity, and position. All LPs have the same total energy E l P = 
ecu, where o) is the statistical weight of the LP (proportional 
to the number of physical particles represented by the LP) 
and € is the energy of each of the physical particles within 
the LP. With all LPs having the same total energy £ LP , the 
largest number of LPs and, thus, the best statistics will be in 
the spectral region where the spectral energy density is 
highest. 

The method of using LPs can be considered to be a gen- 
eralization of the "particle splitting” techniques that are 
routinely used to improve statistics in linear Monte Carlo 
simulations (Pozdnyakov et al. 1983). However, an escape 
probability formalism is not used. Instead, the LPs are split 
as a result of interactions within the system. When an LP 
undergoes interactions, all the physical particles represent- 
ed by the LP are required to experience identical conse- 
quences. Since the total energy of an LP is fixed, € oc 1/w, 
the statistical weight must change as the energy of the 
physical particles, e, changes as a result of the interaction. 
This means that the number of physical particles represent- 
ed by the LP will change. A consequence of this is that each 
interaction requires the elimination or creation of new LPs 
such that energy, momentum, and particle number are con- 
served in a statistical sense. For example, if € doubled 
because of a collision, then the statistical weight of the LP 
would be reduced by a factor of 2. In order to conserve 
particle number, an additional LP, having the same proper- 
ties (e.g., energy, weight, direction) as the first, would be 
created. The reader is referred to Stern et al. (1995a) for 
more details. 

The most important feature of the NLMC method is that 
all the LPs are propagated in parallel. To do this, each LP is 
propagated during each time step, which is chosen to be 
small enough such that <5% of the LPs undergo an inter- 
action. While each LP propagates, the properties of the 
system (e.g., optical depth and temperature) used for deter- 
mining the collision probabilities are determined from the 
state of the system at the end of the previous time step. All 
the properties are updated after the last LP is propagated. 
By using such small time steps, it is ensured that the proper- 
ties of the system will not change significantly during the 
time step, and consequently, each LP will "see” essentially 
the same configuration. Small time steps also allow the 
properties of the system to evolve smoothly. 

For this paper, the types of particles that are simulated 
are photons, electrons, and positrons (although our numeri- 
cal routine also is capable of simulating other particles). 
Although we do not simulate protons, because the cross 
section for Compton scattering is much smaller than that 
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of the electrons and positrons, we specify their initial dis- 
tribution and enforce charge neutrality such that the 
initial distribution of electrons is equal to the specified 
proton distribution. During the simulation, we require that 
n e — n+ = n p , where n e , and n p are the number density 

of electrons, positrons, and protons, respectively. 

The physical interactions that are taken into account are 
Compton scattering, photon-photon pair production, and 
pair annihilation. The numerical code is also capable of 
simulating synchrotron radiation and e~e~ and e~e + 
Coulomb interactions; for simplicity, however, we do not 
consider these processes in this paper. For the reprocessing 
of radiation in the accretion disk, Compton scattering and 
photoabsorption, resulting in fluorescent line emission and 
thermal emission, are considered. Simulations of all “two 
body” processes are carried out using the fully relativistic, 
angle-dependent, quantum electrodynamic cross sections, 
as well as three-dimensional relativistic kinematics 
(Akhiezer & Berestetskii 1965). Continuous processes are 
simulated through discrete “collisions,” using the same 
techniques as for the “two body” processes. However, the 
cross sections that determine the collision rates are chosen 
such that the time-averaged rate of energy losses agrees 
with the continuous, analytic cooling rates. 

In order to allow for nonuniform coronal models, the 
corona is divided into spatial cells. For plane-parallel 
geometry these cells are uniform layers with a thickness A z 
(although it is possible for Az to depend on the cell number). 
Since the radiation field and the properties of the corona 
within each cell are assumed to be uniform, Az is chosen 
such that each cell is optically thin. However, as discussed 
in § 3.3, if the statistical weight within a cell is too small, 
statistical fluctuations become problematic. We find that 10 
spatial cells allow for good resolution of the vertical depen- 
dence of the temperature, opacity, and radiation field while 
also allowing for acceptable statistics within each cell. 

3.2. Modifications to the NLMC Code 

The LPMC code was originally invented to model pair 
cascades in relativistic pair plasmas (Stern 1985). We made 
several modifications to the code to make it more suitable 
for studying accretion disk coronae. These modifications 
include the use of thermal pools and the treatment of radi- 
ation reprocessing in the accretion disk. 

3.2.1. Thermal Pools 

For thermal accretion disk corona models, it is not neces- 
sary to treat all of the electrons and positrons as LPs, since 
it is assumed that most of them are thermalized. Instead, for 
each space cell, a thermal electron “ pool ” and a positron 
“pool” are used. The pools are characterized by having a 
statistical weight w(z), and a temperature T c (z), where i is the 
cell number in the corona. Since bulk motion is not con- 
sidered in this paper, the velocities of the individual 
“particles” within the pool are statistically determined 
from a relativistic Maxwellian distribution appropriate for 
the pool’s temperature. Within each cell, we assume that the 
density distribution of the pool “particles” is uniform. The 
temperature of the positron pool is assumed to be equal to 
the temperature of the electron pool. At the beginning of a 
simulation, all electrons and positrons are assumed to be 
thermalized and are put into their respective pools. When 
an LP is interacting with the pool (e.g., by Compton scat- 
tering with a photon LP), a “miniparticle,” with a weight 


equal to the weight of the interacting LP, is drawn from a 
relativistic Maxwellian distribution. If, as a result of the 
interaction, the “miniparticle” ends up with an energy 
greater than a cutoff energy € max (i), then the “miniparticle” 
is ejected from the pool and converted to an LP. For this 
paper, we set e max (z) = %kT c {i), so that over 99% of the elec- 
trons and positrons are thermalized and only the very ener- 
getic particles are treated as LPs. Alternately, after an 
electron or positron LP interaction or following a pair pro- 
duction event, the electron or positron LP is inserted into 
the corresponding pool if its energy is less than the cutoff 
energy. Once an LP is inserted into the pool, its identity is 
lost, but the pool gains a statistical weight and an energy 
equal to that of the LP. 

The use of thermal pools makes the computations much 
more efficient because it is not necessary to simulate all 
low-energy particle interactions within the pool. These 
physical low-energy interactions thermalize the electrons 
and positrons, and we already have assumed the pools to be 
thermalized. The probability that an LP interacts with the 
pool is determined in the same fashion as for LP-LP inter- 
actions, but the total weight of the pool is used instead of 
the individual weight of the “target” LP. The only pool- 
pool interactions that we simulate are annihilations 
between the positron pool and the electron pool. Pool anni- 
hilations are treated by dividing the positron pool into 100 
“ miniparticles,” each having a weight equal to the total 
weight of the pool divided by 100. Each positron mini- 
particle is given a velocity, drawn from a relativistic Max- 
wellian distribution, and a random direction of motion. 
Each miniparticle is then propagated in time. For each 
annihilation event, a statistical weight wj 100 is taken out of 
both the electron pool and the positron pool, while an 
equivalent weight is inserted in the form of photon LPs. 

3.2.2. Thermal Structure 

The thermal structure of the corona is given by the dis- 
tribution of the pool temperature of each cell, and the non- 
thermal structure is given by the distribution of electron 
and positron LPs. The temperatures of the pools adjust in 
time as a result of energy gains and losses through LP 
interactions, pool annihilations, and external energy 
sources (which we are associating with dissipation of accre- 
tion energy in this paper). The equilibrium pool tem- 
perature is given implicitly by requiring 

tf(z)-C(T,^ + , U r J = 0, (1) 

where H(z) is the external heating rate per unit volume, z is 
the height above the accretion disk, C is the net cooling rate 
per unit volume, n e + are the number densities of electrons 
and positrons, and U Tad is the energy density of the radi- 
ation field. Compton scattering with low-energy photons is 
the dominant method of coronal cooling, and external 
energy deposition is the dominant source of coronal 
heating. 

3.2.3. Disk Reprocessing 

The reprocessing of radiation within the cold accretion 
disk is simulated by using a linear Monte Carlo method. 
Physical processes included are Compton scattering, pho- 
toabsorption, and fluorescence. Since we will apply this 
model to the spectra of BHCs that contain an Fe Ka line, 
we assume that all elements in the disk are neutral and 
simulate Compton scattering off stationary electrons; the 
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main results described below are not influenced by this 
assumption. In this case, the main changes in the reflected 
spectrum are at energies below about 10 keV. We use 
photoionization cross sections from Band et al. (1990) and 
Verner et al. (1993) and solar abundances given by Grevesse 
& Anders (1989) (taking into account the smaller iron abun- 
dance found by Biemont et al. 1991 and Holweger et al. 
1991, where [Fe] = 7.54). Fluorescence in the disk is simu- 
lated by using the fluorescence yields given by Kaastra & 
Mewe (1993). In order to treat multiple scattering and 
absorption/reemission events, we use a variation of the 
method of weights (Pozdnyakov et al. 1983). Photons 
absorbed during a scattering event are reemitted with their 
energy changed to the energy of the fluorescence line and 
their weight reduced by a factor given by the fluorescence 
yield. Energy is conserved by forcing the rate of LPs emitted 
as blackbody photons and line photons to equal the rate of 
absorption events. Photons that scatter within the disk 
more than 50 times are assumed to be thermalized, leading 
to the emission of a blackbody photon. 

3.3. Statistical Fluctuations 

Because of statistical fluctuations that arise from using a 
finite number of LPs in the simulations (currently we use 
N LP = 65,536), the Compton cooling and pair production 
rates vary between iterations, which can cause large tem- 
perature and opacity oscillations. Because of the strong 
coupling between the temperature and the opacity of the 
corona, intrinsic fluctuations in the temperature lead to 
fluctuations in the opacity (and vice versa). These oscil- 
lations can lead to an escaping photon spectrum that differs 
from that of a “steady state” model (even if the two models 
have the same “time averaged” properties). In order to 
reduce the magnitude of the fluctuations, the temperature of 
each electron or positron pool is averaged over the previous 
four iterations. Averaging over a larger number of steps 
leads to overstable oscillations in both the temperature and 
pair opacity, because the corona cannot adjust to changes 
in the radiation field quickly enough, causing it to over- 
shoot its equilibrium values. An additional improvement is 
made by averaging the opacity over the previous 10 iter- 
ations. This averaging still allows for the correct annihi- 
lation and pair production rates, when averaged over all the 
iterations, and therefore does not alter the equilibrium pair 
opacity values. Occasionally, even when averaging, severe 
fluctuations will cause significant temperature changes. We 
prevent this occurrence by artificially preventing the tem- 
perature of the pool to change by more than 10% between 
iterations. When the system is in equilibrium, the fluctua- 
tions of the cooling rates are symmetric around the equi- 
librium values. Thus, even though energy is not conserved 
in the iterations in which large temperature changes are 
artificially prohibited, the energy of the system is conserved 
in a statistical or “ time averaged ” sense. 

An additional source of statistical fluctuations of the pool 
properties is due to unlikely LP-pool interactions. For 
example, there is a very small probability that a Compton- 
scattered photon will give all its energy to a “ pool electron ” 
particle, causing the particle to be ejected from the pool as 
discussed above. For optically thin models, in which the 
statistical weight of the electron pool is small, these ejection 
events can lead to large fluctuations in the opacity. In order 
to reduce the magnitude of these fluctuations, LPs that 
interact with a pool are “ split ” into N m mini-LPs, where the 


mini-LPs have a statistical weight equal to \/N m of the 
original weight of the LP. Each of the mini-LPs is then 
propagated in the usual way, and the final state of the LP 
and pools is determined in a statistical sense. We find that 
N m = 25 allows for adequate statistics without degrading 
the efficiency of the code too much. For models with a total 
coronal optical depth oft < 0.1, we let N m = 50. 

Using the above methods, the rms statistical fluctuations 
of both the coronal temperature and the total optical depth 
are a rms < 5% after the system has reached its “steady 
state” equilibrium values. We emphasize that, although the 
simulations are integrated in “time,” the true time- 
dependent behavior of the simulations may not be correct, 
since the statistical methods force the timescales for the 
changes in the temperature and opacity to be longer than 
their proper values. Thus, time-dependent simulations are 
only meaningful on timescales longer than the time needed 
for the coronal properties to vary smoothly. In this paper, 
we study only “steady state” models and use “time” only 
insofar as it measures the approach to an equilibrium. Thus 
the transient behavior (which arises from simulating one 
model by starting from the ending point of a different 
model) of the radiation field and the coronal properties are 
disregarded. Only after the simulation has reached a steady 
state do we begin recording the spectrum and determining 
the corona properties. 

4. SLAB-GEOMETRY ADC 

Although more complicated geometries are easily simu- 
lated with our code, for this paper we consider only the 
plane-parallel geometry. In addition, although we can have 
nonuniform proton density distributions and nonuniform 
heating rates, here we consider uniform models so that we 
can compare our results with previous work. Last, we 
assume that the protons have the same temperature as that 
of the electrons. In future papers, we will relax all of these 
simplifications. 

In our models, the rest-mass energy density of electrons is 
several orders of magnitude lower than the radiation energy 
density. Therefore the collision rate between photons and 
electrons is much higher than the Coulomb collision rate. 
The higher Compton collision rate indicates that the radi- 
ation field and the cooling rate of the corona are dominated 
by Compton collisions, so we neglect cooling due to 
thermal bremsstrahlung. For simplicity, we also neglect 
synchrotron radiation for this paper. 

In the plane-parallel (slab) geometry, the accretion disk 
corona is above and below the optically thick (but geo- 
metrically thin) accretion disk, as shown in Figure 1. For 
this model, the geometric covering factor Q/2n = 1. Similar 
to HM93, we assume that all corona and disk properties are 
constant with respect to radius and assume azimuthal sym- 
metry. The principal parameters of our models are (1) the 
seed optical depth of the corona, t p (which does not include 
the contribution due to pairs), (2) the local compactness 
parameter of the corona, / c , and (3) the temperature of the 
accretion disk, T BB . For slab geometry, we define a local 
compactness parameter of the corona as 


where z 0 is the scale height of the corona, rn e is the mass of 
the electron, c is the speed of light, and is the rate of 
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Fig. 1. — Plane-parallel (slab) geometry. The accretion disk is assumed 
to be optically thick; all photons that are incident onto the disk are either 
reflected or photoabsorbed. The disk is assumed to be symmetric around 
its midplane ( hatched region). 


energy dissipation per unit area into the corona. We remind 
the reader that l c is a local compactness parameter and is 
related to the global compactness parameter used by HM93 
by /(local) = / c (global)/7T for cylindrical disk geometry with 
a uniform 4/.. The optical depth of the corona is given by 
t t = jo Gjn(z)dz , where <r x is the Thomson cross section 
and n{z) is the number density of all electrons and positrons. 

4.1. Energetics 

We define a local disk compactness parameter as 


1 r- 

d ~ m„ c 3 Z ° 


(3) 


where F d is the flux of radiation emitted by the disk (ergs 
s“ 1 cm -2 ) in the form of thermal radiation and fluorescence 
lines. Radiation that is reflected off the disk is not included 
in F d . The flux of radiation that is absorbed by the accretion 
disk is given by 

F abs = F c -(l -<i), (4) 

where F~ is the flux of downward-directed radiation within 
the corona just above the accretion disk and a is the albedo 
of the disk (averaged over all energies and angles). We 
require that the flux of radiation emitted by the disk equal 
the sum of the rate of gravitational energy dissipated within 
the disk and the flux of downward-directed radiation 
absorbed by the disk. This requirement is expressed as 


F d = (1 -f)P G + F a 


(5) 


where P G is the rate of gravitational energy dissipated per 
unit area and /is the fraction of this energy that is deposited 
directly within the corona. By defining the total com- 
pactness parameter due to gravitational dissipation to be 
Iq = (<r T /m e c 3 )z 0 P G , equation (5) becomes 


i* = « * -/) + l 


( 6 ) 


where / abs = (a T /m e c 3 )z 0 F abs . Finally, by defining the 
escaping radiation compactness parameter as / esc == 
(<7 T /m e c 3 )z 0 F esc , the requirement of conservation of energy 
can be expressed as 

= Ksc = K + L ~ Lbs - ( 7 ) 

Unlike HM93, we do not assume /= 1. Instead, we allow 
for an intrinsic disk compactness by setting (1 — f)l G = 1 
for all models. The fraction of gravitational energy dissi- 


pated into the corona is then given by 

L 


f= 


1 + /< 


( 8 ) 


If l c > 1, then / ~ 1 since most of the gravitational energy is 
dissipated into the corona and F d % F th . Conversely, if l c < 
1, then f zz 0 since most of the gravitational energy is dissi- 
pated directly in the disk and F d « P G . We study models for 
0.01 < l e < 1000, corresponding to 0.01 </< 1. Although 
setting l G ( 1 — /) to unity is an arbitrary choice, the main 
results presented in this paper are independent of the value. 
For a given total optical depth, the self-consistent coronal 
temperature depends on the ratio of the coronal com- 
pactness parameter to the intrinsic disk compactness 
parameter, which is a function of / (for a more detailed 
discussion, see HM93). 

The temperature of the disk, ?BB> is held fixed, as in 
HM91. The spectral shape of the thermal radiation is given 
by Planck’s law, normalized such that the total flux of 
thermal radiation is consistent with the self-consistent com- 
pactness parameter of the disk (including reprocessing). 
With a fixed T BB , different values of l d correspond to differ- 
ent values of the coronal scale height z 0 . 

With / c , r e , and 7^ B specified, the NLMC code is used to 
find the resulting temperature structure, T(z), the total 
opacity, t t (including the contribution by pairs), the internal 
radiation field, and the escaping radiation field. The spec- 
trum of escaping radiation is stored in 10 bins in g between 
0 and 1, where g — cos 6 and 6 is the angle between the 
normal of the disk and the photon’s direction. We also 
compute the spectrum of radiation incident onto and reflec- 
ted from the disk, also stored in 10 bins in g. As mentioned 
in § 2, the reprocessing of radiation in the disk is computed 
by interfacing a linear MC technique with our LPMC code. 
We are able to calculate the energy-integrated disk albedo, 
a, more accurately than HM93 since we do not assume that 
the radiation field incident onto the disk field is isotropic; 
the albedo does depend on the angular distribution of the 
incident radiation field. We stress that we do not simply add 
the reflected component of the radiation, modified by an 
escape probability, to the escaping spectrum. Instead, we 
propagate the reflected component through the corona in 
the same fashion as all other components of the radiation 
field. Therefore, multiple reflections and interactions of 
reflected photons with the corona are taken into account. 

4.2. Generation of the Models 

We have produced a grid of roughly 100 ADC models 
with a slab geometry. For all cases, the spatial distributions 
of the energy dissipation rate, H(z\ and the seed opacity 
distribution, n p (z\ are uniform. For the purposes of compar- 
ing our simulations with the spectra of BHCs (which we do 
in a companion paper, Dove et al. 1997, hereafter Paper II), 
we allowed the seed opacity to vary within the range 0.05 < 
z p < 2.0 and the compactness parameter to vary within the 
range 0.1 <l c < 10 3 . Except for the soft X-ray excess, the 
output spectra are not sensitive to the value of and 
therefore we have fixed kT m = 200 eV (see discussion in § 5). 

The grid is produced by sweeping the coronal com- 
pactness parameter from its minimum to its maximum 
values while keeping the seed opacity fixed. After each 
sweep, the seed opacity is incremented. Each new simula- 
tion is started using the end state of the previous simulation. 
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In this approach, CPU time is saved since the system does 
not have to evolve much prior to reaching an equilibrium 
state as compared with starting a simulation from scratch. 
The only caveat regarding production of the grid is that the 
relative increments in l c should not be more than about 
threefold; for larger increments, the system is unable to 
smoothly evolve to the new equilibrium solution since the 
transitional behavior will heavily disturb the system. We 
have verified that the results presented here are independent 
of the order in which the grid is produced (no hysteresis). 
After the system has reached an equilibrium, the spectrum 
of escaping radiation is binned and the self-consistent 
properties of the corona are determined. 

5. RESULTS 

5. 1 . Consistency Tests and Comparisons with Linear 

MC Models 

In order to test our model, we compared the spectra of 
escaping radiation, produced by our linear version of the 
model (where we do not include reprocessing, pair pro- 
duction, or annihilation, and we use uniform density and 
temperature structures), with the linear MC models of 
Wilms et al. (1996) and J. G. Skibo (1996, private 
communication). Both slab and spherical geometries were 
considered. In all cases, our models are in excellent agree- 
ment with all of the linear MC models. In addition, we 
computed the amplification factor A , where A — L om /L input , 
with L out and L input being the luminosity (ergs s" 1 ) of the 
escaping radiation and the luminosity of the injected seed 
photons, respectively. We have compared our amplification 
factors with those given by Gorecki & Wilczewski (1984) for 
spherical geometry. These amplification factors agree to 
within 7%, and we attribute the differences to numerical 
fluctuations. Because of computational speed constraints, 
Gorecki & Wilczewski did not simulate many particles. 
Amplification factors were also computed using the linear 
MC code of Wilms et al. (1996). Our values agreed with 
these values to within 2%. As an additional test, we com- 
puted the average number of scatterings that a photon 
undergoes prior to escaping the system. We compared these 
results with the analytic expectations for both the optically 
thin and optically thick cases. We also compared the 
average change in the photon’s energy per scattering event, 
as well as the average change in the square of the energy 
change «e 2 », and compared them with the analytic expec- 
tations for both the relativistic and the nonrelativistic cases. 
In all cases, our results are perfectly consistent. 

We compared the spectrum of reprocessed (i.e., reflected) 
radiation when injecting the cold accretion disk with a 
power law with the spectra given by George & Fabian 
(1991) and Haardt (1993). Finally, we compared our 
maximum self-consistent coronal temperature values, as a 
function of opacity, with those given by Stern et al. (1995b), 
Poutanen & Svensson (1996), and Haardt et al. (1997). 
Although our results agree with Stern et al. (1995b) for all 
optical depths, there are some discrepancies with the other 
models for t t < 0.2. For example, for t t = 0.2, we predict 
an average coronal temperature T c = 138 keV while Haardt 
et al. (1997) predict T c = 160 keV, a discrepancy of 12%. 
The discrepancy increases as the opacity decreases. We do 
not understand these discrepancies, and an investigation 
that will address this problem is currently underway. Since 
the discrepancy is less than 10% for x T > 0.2, however, these 


differences are not important for understanding the hard 
spectra of BHCs, as all “best fit” slab models have a total 
opacity t x > 0.3. 

5.2. Self-consistent Thermal Properties of the Corona 

In Figure 2, we show the self-consistent temperature as a 
function of the total optical depth, for several values of l c 
and seed opacities. The solid line marks the maximum 
coronal temperature, as a function of the total opacity, 
T mAX (x T ). This figure is similar to Figure 1 a of HM93 and 
Figure 1 of Stern et al. (1995a), but for clarity we include 
models that are both pair dominated and non pair domi- 
nated, showing how the temperature evolves with increas- 
ing compactness parameter. As pointed out by HM93, the 
maximum self-consistent temperature, for a given total 
opacity, is independent of the seed opacity. However, for a 
given total opacity, pair-dominated models reach the 
maximum temperature with a coronal compactness param- 
eter that is much lower than the corresponding value for 
non-pair-dominated models. 

Also in Figure 2 (inset), we show how the temperature of 
the corona varies with the coronal compactness parameter, 
starting with a seed opacity of z p = 0.2. Note that, for small 
l c values, the temperature increases as a function of l c while 
the opacity remains nearly constant (pair production is 
negligible). Here the Compton cooling rate is dominated by 
the soft photons implicitly emitted by the cold accretion 
disk (/<^ 1), and an increase in l c corresponds to an increase 
in the heating rate without much of an increase in the 
cooling rate. As the corona becomes hotter, the pair pro- 
duction rate increases. Once pair production begins to be 
significant, the increase in t t results in an increase of radi- 
ation reprocessing within the cold accretion disk. Since 
most ( ~90%) of this reprocessed radiation is reemitted as a 
thermal blackbody, the increase in the reprocessed radi- 
ation causes a very large increase of the Compton cooling 
rate within the corona. Thus, through the production of 
pairs and the increase in reprocessed radiation, the 
Compton cooling rate increases more than linearly with 
increasing l c . Since the heating rate is only proportional to 
/ c , the coronal temperature reaches a maximum value and 
then decreases with increasing l c for models in which the 
seed electron optical depth is held constant. We note that 
this T c (l c ) relationship is due to our assumption that the 
intrinsic compactness parameter of the disk is constant 
while l c varies. However, this behavior does yield physical 
insight about how the onset of pair production forces the 
temperature to reach a maximum value and then decrease 
with increasing values of the coronal compactness param- 
eter. Any model in which the flux of the disk becomes domi- 
nated by reprocessing of coronal radiation will exhibit this 
behavior. 

For the models presented here, the pair-dominated 
models do not predict an observable annihilation line in the 
escaping spectrum; thus, for a given total opacity and 
coronal temperature, the pair-dominated and non pair- 
dominated models predict the same spectrum of escaping 
radiation. The degeneracy of models is shown in Figure 2, 
where models with different seed opacities and compactness 
parameters sometimes have the same total optical depth 
and temperature. For this reason, our results are indepen- 
dent of the value of the intrinsic disk compactness param- 
eter, / c (l — /), which we set to unity. If this value were set to 
a higher value, the Compton cooling rate would be higher, 
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Total Optical Depth 

Fig. 2. — Allowed temperature and opacity regime for self-consistent ADC models with a slab geometry. The solid line is derived from a “ fit by eye ” to the 
numerical results. For a given total optical depth, temperatures above the solid line are not possible. The hatched region below this line marks the parameter 
space in which the contribution to the total opacity by pairs is significant. For all models, the blackbody temperature of the disk is kT m = 200 eV. Different 
symbols represent models with different seed opacities. Inset : Average coronal temperature and total opacity for several values oil c . 


but the same maximum temperatures would be reached 
with higher l c values. If / G (l — /) were decreased, the models 
in which l c 1 would not be altered, since the thermal 
emission is already dominated by reprocessed radiation. 
For models in which l c < 1, decreasing / G (l — /) while 
holding the seed opacity constant simply reduces the 
Compton cooling rate and precludes the cooler regions of 
Figure 2 (lower left) from being self-consistent (i.e., reducing 
the value of the intrinsic disk compactness parameter drives 
the solutions toward the pair-dominated, /~ 1, regime). All 
of these arguments have been confirmed by numerical simu- 
lations. 

The blackbody temperature used for our models is T BB = 
200 eV, a value believed to be appropriate for BHCs (see 
Paper II). However, the results presented in Figure 2 are 
insensitive to the disk temperature. For example, a corona 
with a total optical depth r T = 0.25 and a disk temperature 
T bb = 5 eV has a maximum temperature T max =117 keV, 
compared to T max = 128 keV for the T m = 200 eV case. 
These values are consistent with the results of HM93. 
Therefore, models with very low disk temperatures, which 
are attractive in the sense that the thermal excess might not 
be observable as a consequence of the efficient Galactic 
absorption of ultraviolet radiation, have maximum coronal 
temperatures that are colder than the values presented 
above. On the other hand, models with a higher blackbody 
temperature predict even larger thermal excesses than those 
corresponding to a 200 eV disk. 

We conclude that accretion disk coronae cannot have a 
temperature higher than ~140 keV if the total optical 
depth is >0.2, regardless of the coronal compactness 
parameter, intrinsic accretion disk compactness parameter, 
and blackbody temperature. Alternatively, ADCs with a 
temperature >150 keV cannot have an optical depth 


greater than ~0.15 (see Fig. 2). The range of allowed self- 
consistent values of the coronal temperatures and opacities 
for the slab geometry corresponds to a minimum value of 
the angle-averaged photon index T — 1.8 (as determined by 
numerically fitting the spectra between 5 and 30 keV with a 
power law), which places severe limitations on the applica- 
bility of these ADC models to BHCs. However, because of 
the anisotropy break (HM93), the power-law portion of the 
spectrum hardens as the inclination angle decreases (see 
Fig. 5 below). For models with r T = 0.3 and kT c = 1 10 keV, 
where the angle-averaged index T « 1.8, T varies between 
2.1 and 1.7 as the orientation varies from “edge on” to 
“ face on,” respectively. 

In the accompanying Paper II, we demonstrate that ADC 
models with a slab geometry cannot explain the broadband 
X-ray spectrum of Cyg X-l and other BHCs. The observed 
hard power laws in BHCs can be explained only with 
models that have photon-starved sources, where Compton 
cooling is less efficient and the coronae are allowed to have 
higher temperatures. We find that ADCs with a spherical 
geometry and an exterior cold accretion disk allow for such 
conditions. 

5.3. Nonuniform Temperature Structure 

As discussed in § 4, our models allow for the determi- 
nation of the temperature and opacity within local cells, 
which allows for the possibility of nonuniform temperature 
structures. In Figure 3, we give the relative temperature 
deviation (A T/<T C » for assorted values of t p , T c , and / c . 
Note that, for optically thin models (t t < 0.3), the tem- 
perature is essentially uniform, while the temperature varies 
by up to a factor of 3 for optically thick models. For t t > 
0.5, the region nearest the accretion disk is always the 
coldest, while the region farthest away from the disk (largest 
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Fig. 3. — Vertical temperature structure of the accretion disk corona 
expressed as the relative deviation of the local temperature, T(z), from the 
volume-averaged coronal temperature (T c >. Several models, with different 
combinations of the seed optical depth and coronal compactness param- 
eter, are shown. 


values of z) is the hottest. This is easily understood since the 
energy density of the radiation field decreases with increas- 
ing height, resulting in a decrease of the Compton cooling 
rate with z. For optically thin coronae, the energy density of 
the radiation field does not vary too much with height; 
however, the average photon energy of the radiation field 
increases with height since the spectrum of the internal radi- 
ation field hardens with increasing height. Therefore, since 
the Compton cooling rate is proportional to both the 
average energy, <€>, and the energy density, the minimum 
temperature can occur at intermediate heights rather than 
always at z = 0. 


5.4. The Spectrum of Escaping Radiation 
5.4.1. Comparison with Uniform Models 

In Figure 4, for t t = 0.25, 0.50, and 2.0, we compare the 
escaping spectra between models with a uniform tem- 
perature structure and nonuniform models (arising from a 
uniform heating distribution). For the optically thin models, 
the thermal gradients do not give rise to any significant 
changes in the spectral shape of the radiation field. In fact, 
as shown in Figure 4, the spectra are virtually indistinguish- 
able, as the fractional difference between the uniform and 
nonuniform models is no more than 1% for the optically 
thin models. Therefore, for optically thin ADC models, a 
nonuniform distribution of the heating rate is required to 
sustain temperature gradients large enough such that the 
spectrum of escaping radiation significantly differs from the 
uniform models. For parameters appropriate for low-mass 
X-ray binaries, where t T > . 1, the spectrum of the escaping 
radiation is significantly different from the corresponding 
spectrum resulting from a uniform model. The nonuniform 
model produces a slightly harder high-energy tail and is 
caused by the small contribution of the uppermost region of 
the corona, where the temperature is hotter than the 
average coronal temperature. 

5.4.2. Anisotropy of the Escaping Radiation Field 

In Figure 5, we show the spectrum of the escaping radi- 
ation for several values of the inclination angle. We define 
p = cos 0, where 9 is the angle between the line of sight and 
the normal of the accretion disk. As expected, the spectral 
shape does not vary with respect to the inclination angle for 
optically thick models. However, for optically thin models, 
the spectral shape of the escaping radiation field (and the 
internal radiation field) does vary with p. Since the effective 
optical depth for radiation leaving the accretion disk is x jp. 
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Fig. 4. — Comparison of the predicted spectra of escaping radiation from our nonlinear ADC models with that of our linear models, where the vertical 
temperature distribution is uniform. From top to bottom, the average temperature and total optical depth (k(T c } [keV], t t ) are (11 8.4, 0.25), (69.8, 0.50), and 
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Fig. 5.— Spectra of escaping radiation for several values of /x, where 
fi = cos 0. The uppermost plot is the spectrum averaged over all angles. 
Below, in order, fi = 0.9, 0.7, 0.5, 0.3, and 0.1. For all cases, the average 
coronal temperature is k(T c y = 118 keV, the optical depth is t t — 0.25, 
and the blackbody temperature of the disk is kT BB = 200 eV. 


the fraction of blackbody radiation that escapes the system 
without subsequent interactions with the corona decreases 
with decreasing \i. Consequently, the thermal excess 
becomes weaker with decreasing values of p (these argu- 
ments also apply to the strength of the Fe Ka fluorescence 
line). Finally, the spectrum becomes softer as p decreases 
since the magnitude of the reflection “hump” decreases 
with the effective optical depth. These results are in agree- 
ment with HM93. Therefore the spectra of ADCs that are 
seen “face on” (^ — 1) are the hardest, but they also have 
the largest amount of thermal excess. On the other hand, 
ADCs viewed “edge on” do not have an observationally 
recognizable soft excess in their spectra, nor do they contain 
any reprocessing features. If the matter responsible for 
producing the “reflection features” is the optically thick, 
cold accretion disk that is also producing the seed photons, 
it does not appear possible that the spectra of ADCs can 
have strong reflection features without also containing a 
strong soft excess. While a very cold accretion disk (i.e., 
/cT bb < 10 eV) would produce thermal emission that could 
be “ hidden ” because of the efficient Galactic absorption of 
ultraviolet radiation, it is questionable whether an accretion 
disk of a BHC, in the radial region where most of the seed 
photons are produced, could be so cold (for a more detailed 
discussion, see Paper II). 

5.5. Strength of Fe Ka Fluorescence Line 

The strength of the Fe Ka fluorescence line for BHCs 
provides stringent constraints on the properties of ADC 
models. For a slab geometry, roughly half the Comptonized 
radiation field within the corona is reprocessed within the 
cold accretion disk. For parameters appropriate for BHCs, 
where the radiation field is hard (having a photon index of 
T ~ 1.5), the large amount of reprocessing gives rise to a 
very large equivalent width (EW) of the Fe Ka line. For 
optically thin models having a temperature 130 keV > 
kT c > 100 keV, the angle-averaged EW of the iron line is 
120 eV > EW > 90 eV (for each temperature, there is a 
scatter of about 10 eV due to the range of optical depths, the 
uncertainty in the EW measurements, and the statistical 
nature of the Monte Carlo simulations). The predicted EWs 
would be even higher than these values if a harder spectrum 


(i.e., a spectral shape that describes the Cyg X-l obser- 
vations better than the self-consistent spectrum) were inci- 
dent onto the cold accretion disk. For example, according 
to simulations with our linear Monte Carlo code, a radi- 
ation field with a power-law index of 1.5 and an exponential 
cutoff at 150 keV irradiated onto the cold disk will result in 
an EW — 150 eV for a “face on” orientation. Although 
there is a great deal of uncertainty, the observed EW of the 
Fe Ka line of BHCs is very small, usually having only an 
upper limit rather than a true detection. For Cyg X-l, the 
EW is <70 eV (Ebisawa et al. 1996; Gierlinski et al. 1997). 
Therefore, unless the abundance of iron is much less than 
solar, the slab-geometry models predict EW values that are 
significantly larger than the values observed for Cyg X-l. As 
we discuss in Paper II, this problem is not encountered with 
models having a sphere-plus-disk geometry. For that 
geometry, the accretion disk, as seen by the corona, has a 
covering factor roughly one-third of that for the slab 
geometry, giving rise to less reprocessing of coronal radi- 
ation (Done et al. 1992). 

5.6. Comparisons with Previous Models 
5.6.1. Hua & Titarchuk {1995) 

The major difference between our models and the analy- 
tic models of Sunyaev & Titarchuk (1980) and Hua & Titar- 
chuk (1995) is that our models take into account the 
reprocessing of radiation within the cold accretion disk, 
which gives rise to reflection features in the spectra. In addi- 
tion, our models give a more accurate description of the 
transition between the thermal blackbody (seed photons) 
and the Comptonized portion of the spectrum, since Titar- 
chuk’s analytic model is only valid for energies much higher 
than the energy of the seed photons (nor does this model 
include any anisotropic effects). Therefore, it is expected 
that they will predict different spectra of escaping radiation. 
However, because of the popularity of this model, we give a 
comparison (Fig. 6) of our models with the XSPEC version 
of Titarchuk’s slab-geometry Comptonization model (Hua 
& Titarchuk 1995). Here the spectral shape of the seed 
photons for Titarchuk’s model is described by a Wien law 
rather than a Planckian distribution. Although the two 
models differ significantly, they do agree for energies much 
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Fig. 6. — Comparison of the predicted spectrum of escaping radiation 
with the analytic model of Titarchuk (Hua & Titarchuk 1995), for several 
inclination angles. Solid lines correspond to our model, and the dashed 
lines correspond to Titarchuk’s model. From top to bottom, fi = 0.1, 0.5, 
and 0.9. Again, k<T c > = 1 18 keV, t T = 0.25, and kT BB = 200 eV. 



No. 2, 1997 


THERMAL ACCRETION DISK CORONA MODELS. I. 


757 



10 100 1000 
Energy [keV] 

Fig. la 


CO 

CC 


1.4 


1.2 


1.0 

'f 

. j 1 

0.8 


j . 


m 


r p 


if 


10 100 
Energy (keV) 

Fig. lb 



Fig. 7. — Comparison of the predicted spectrum of escaping radiation with spectra of the form F(E) oc E exp( EfE cut ). (fl) The uppermost spectrum is 
the angle-averaged spectrum; then, from top to bottom, p. = 0.1, 0.3, 0.5, 0.7, and 0.9. The model parameters are the same as in Fig. 6. (b) Residuals of the 
numerical model divided by the analytic model for p = 0.9. 


higher than kT BB if the reflection features are ignored. 
However, our high-energy tails do not decrease with 
increasing energy as rapidly as the analytic model. This 
disagreement has been discussed by several other authors 
(Stern et al. 1995a; Skibo et al. 1995). We also note that, 
with the implementation of XSPEC available at the time of 
this writing (version 10.0), it is not possible to add a reflec- 
tion component to Titarchuk’s Comptonization model self- 
consistently. 

5.6.2. Power Law with Exponential Cutoff 
As shown in Figure la , we have attempted to describe 
our simulated spectra of escaping radiation by a power law 
with an exponential cutoff, 

N e = F 0 E~ r exp ( — £/£ 0 ) , (9) 

where N E is the photon flux as a function of energy, E. 
Although we were able to achieve good “ by eye ” fits for the 
angle-averaged spectrum (integrated over all values of p), 
we were unable to find good fits for individual bins of p. 
This problem is due to the presence of reflection features. In 
Figure lb , we show the residuals of our numerical model 
divided by the analytic model (eq. [9]) for p = 0.1. The fits 
become better as the inclination angle increases, since the 
strength of reprocessing features decreases with increasing 
angle (§ 5.4.2). It appears to be a coincidence that, for the 
angle-averaged spectrum, the reflection features add up in 
such a way that a power law is better preserved. 

5.6.3. Reflected , Exponentially Cut Off Power Law 
Using XSPEC, we have compared our model with the 
reflected, exponentially cut off power-law model PEXRAV 
(Magdziarz & Zdziarski 1995). Using the Rossi X-Ray 
Timing Explorer ( RXTE ) proportional counter array 
(PCA) response matrix, we simulated a 10 4 s observation 
using FAKE1T and our NLMC model, where kT BB = 200 
eV, k(T c } =118 keV, t x = 0.25, and p = 0.9. As shown in 
Figure 8, we fitted these “data” with a superposition of a 
200 eV blackbody, a Gaussian line at 6.4 keV, and 
PEXRAV. Our best fit yielded a photon index of 2.0, a 
cutoff energy E c = 400 keV, and the relative reflection 
parameter/ = 0.48. We froze the abundance parameter to 
unity and set p = 0.9. It is clear that there is much more 
structure in our NLMC model than in the reflection model. 


showing a pitfall of assuming that a Comptonized, unre- 
flected spectrum can be described by an exponentially cut 
off power law. It is also interesting that the best-fit covering 
fraction is/< 0.5, given that the physical covering fraction 
for the slab geometry is/ = 1. 

6. CONCLUSIONS 

We present accretion disk corona models in which the 
radiation field, the temperature, and the total opacity of the 
corona are determined self-consistently. We take into 
account the coupling between the corona and the accretion 
disk by including reprocessing of radiation in the accretion 
disk. 
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Fig. 8.— Comparison of the reflected, exponentially cut off power-law 
model, PEXRAV (Magdziarz & Zdziarski 1995), with the NLMC model. 
Using XSPEC (Arnaud 1996), “fake” data were produced by folding the 
NLMC model through the RXTE PCA response matrix. The parameters 
of the NLMC model are the same as in Fig. 6. A Gaussian at 6.4 keV and a 
blackbody distribution with JcT BB = 200 eV have been added to the 
PEXRAV model, in which an isotropic radiation field given by 
N(E) oc E~ t exp {-E/E c ) is reflected by a slab of cold matter. The “best 
fit” reflection model resulted in an equivalent width of the Fe Ka line 
EW « 100 eV, T % 2.0, E c - 400 keV, and the covering fraction of the disk 
/= 0.48. Galactic absorption has been taken into account by both models 
using WABS (using the Morrison & McCammon 1983 cross sections) and 
assuming N H = 6 x 10 21 cm 2 . 
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The range of self-consistent temperatures and total opa- 
cities places severe limitations on the model’s applicability 
in explaining the X-ray spectra of BHCs. The maximum 
coronal temperature for models having t x > 0.3 is kT c % 
110 keV. The corresponding angle-averaged, 5-30 keV 
photon index is T > 1.8, a value too large to explain the 
observed hard spectra of BHCs (with typical photon power- 
law indices of ~ 1.5). Even models with a “ face on ” orienta- 
tion predict photon indices of T > 1.7. The spectra (for 
energies E > 5 keV) of BHCs, most notably Cyg X-l, have 
been adequately described by linear ADC models in which 
r ~ 0.3-0.5 and kT ~ 150 keV (e.g., Haardt et al. 1993; 
Titarchuk 1994). None of these previous ADC models are 
self-consistent, for they use temperature and opacity values 
that are outside the allowed region. 

In addition, as described in Paper II, attempts to describe 
the broadband spectrum of Cyg X-l have proved to be 
much more difficult than modeling the spectrum over a 
smaller energy range. This difficulty is due to the fact that 
optically thin ADC models always predict a very large 
thermal excess, a feature not found in the observations, 
unless the inclination angle is nearly “edge on.” For the 
case of Cyg X-l, these high inclination angles probably can 
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spectra. Optically thick models, which also predict a small 
thermal excess, have a self-consistent coronal temperature 
that is much too cold to explain the spectra of BHCs. As a 
consequence of Galactic absorption, the predicted thermal 
excess decreases with decreasing values of the blackbody 
temperature. However, a colder accretion disk results in a 
higher Compton cooling rate, yielding even lower coronal 
temperatures than the values corresponding to a 200 eV 
disk. 
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